/*
 * mplib.c
 *
 *  Created on: Apr 7, 2010
 *      Author: nguyencao
 */


#include "mplib.h"
#include "utils.h"
#include "math.h"

/**
 * Read a line from the stream of file
 * return 0 if no error
 * 		  else if EOF or error return ERROR
 */
int getlin(FILE *file, CHAR **line) {
	if(file == 0) {
		return ERROR;
	}
	if(fgets(*line, MAXLINE, file) == NULL) {
		if(feof(file))
			return EOF;
		else
			return ERROR;
	}
	return 0;
}

/**
 * Put a line to the file
 * Return 0 if success
 * 		  else return ERROR if not
 */
int putlin(CHAR *line, FILE *file) {
	if(file == 0) {
		return ERROR;
	}
	if(fputs(line, file) == EOF) {
		return ERROR;
	}
	return 0;
}

/**
 * Output error on the standard output
 */
void error(CHAR str[]) {
	printf("\n%s",str);
	exit(1);
}

/**
 * Output a large integer n into file
 * Break n into subgroup, each with maximum length is w
 */
void putdec(FILE *file, INTEGER n, int w) {
	CHAR *ch = 0;
	int i,j;
	int nd = itoc(n, &ch, 10, 79);
	for(i=nd;i<w;i++)
		j=fputc(BLANK, file);
	for(i=0;i<nd;i++)
		j=fputc(ch[i],file);
}

/**
 * Output a small integer into file with maximum length is w
 */
void putdec2(FILE *file, int n, int w) {
	CHAR *ch = 0;
	int nd = itoc2(n,&ch, 10, 79);
	int i,j;
	for(i=nd;i<w;i++)
		j=fputc(BLANK, file);
	for(i=0;i<nd;i++)
		j=fputc(ch[i],file);
}

/**
 * Convert large integer n to char string in str
 * Return the size of the str if < size
 * or size otherwise
 */
int itoc(INTEGER n, CHAR** str, int base, int size) {
	*str = mpz_get_str((CHAR *) 0, base, n);
	if(strlen(*str) > size) {
		return size;
	}
	return (strlen(*str));
}

/**
 * Convert small integer n to char string in str
 * Return the size of the str if < size
 * or size otherwise
 */
int itoc2(int n, CHAR** str, int base, int size) {
	INTEGER N;
	init2(N,n);
	*str = mpz_get_str((CHAR *) 0, base, N);
	if(strlen(*str) > size) {
		return size;
	}
	return (strlen(*str));
}

/**
 * Convert the string starting at position i of str
 * to large integer, then increment i
 */
void ctoi(CHAR *str, int *i, INTEGER k) {
	int j = *i;
	int start, end;
	CHAR *s;

	while(str[j] == BLANK)
		j++;
	start = j;
	while(isdigit(str[j]))
		j++;
	end = j;

	s = (CHAR *)malloc(sizeof(CHAR)*(end-start + 1));
	s[end - start] = EOS;
	for(j=start;j<end;j++) {
		s[j-start] = str[j];
	}
	mpz_set_str(k,s,10);
	*i = end;
}

/**
 * Convert the string starting at position i of str
 * to small integer, then increment i
 */
int ctoi2(CHAR *str, int *i) {
	int j=0;
	int d = 0;
//	printf("A2.1.1.\n");fflush(STDOUT);
	while(str[*i] == BLANK)
		(*i)++;
//	printf("A2.1.2 i=%d\n",*i);fflush(STDOUT);
	for(j=0;str[*i] != EOS; *i=*i+1) {
		d = (int)str[*i];
		if(d<48 || d>57) {	// 48: digit 0, 57: digit 9
			break;
		} else {
			d=d-48;
		}
		j = 10*j+d;
	}
	return j;
}

// copy A into B: B = A
void mpcopy(INTEGER A, INTEGER B) {
	mcopy(B,A);
}

void mult1(INTEGER a, int b, INTEGER c) {
	char *negb = " negative b in mult1 ";
//	mpchk(400,a);
	if(b > 0) {
		INTEGER b2;
		init2(b2,b);
		mul(a,b2,c);
//		mpchk(402,c);
		clear(b2);
	} else if(b < 0) {
		putdec2(STDOUT, 401, 6);
		putdec2(STDOUT, b, 10);
		putlin(negb, STDOUT);
		putc(NEWLINE, STDOUT);
	}
	set2(c,0);
}

int div1(INTEGER a, int b, INTEGER q) {
	char *negb = " b <= 0 in div1 ";
	int r = 0;
//	mpchk(400,a);
	if(b > 0) {
//		mpchk(500,a);
		div_qr(a,b,q,&r);
	} else if(b < 0) {
		putdec2(STDOUT, 501, 6);
		putdec2(STDOUT, b, 10);
		putlin(negb, STDOUT);
		putc(NEWLINE, STDOUT);
		exit(1);
	}
	return r;
}

// check dimension and digits of the large integer a
void mpchk(int code, INTEGER a) {
	char *s5 = " mpchk ";
	char *s1 = " dimension < 2 ";
	char *s2 = " dimension > MPDIM ";
	char *s3 = " negative digit ";
	char *s4 = " digit >= MPBASE ";
	int *arr = 0;
	int d = 0, i=0;
	convert(a,&arr,&d);
	putlin(s5,STDOUT);
	putdec2(STDOUT,code,5);
	for(i=0;i<d;i++) {
		putdec2(STDOUT,arr[i],10);
	}
	putc(NEWLINE, STDOUT);
	if(d < 2) {
		// dimension is too small
		putdec2(STDOUT,code,5);
		putdec2(STDOUT,d,10);
		putlin(s1,STDOUT);
		putc(NEWLINE, STDOUT);
		exit(1);
	} else if(d > MPDIM) {
		putdec2(STDOUT,code,5);
		putdec2(STDOUT,d,10);
		putlin(s2,STDOUT);
		putc(NEWLINE, STDOUT);
		exit(1);
	} else {
		for(i=0;i<d;i++) {
			if(arr[i] < 0) {
				// negative digit
				putdec2(STDOUT,code,5);
				putdec2(STDOUT,d,10);
				putdec2(STDOUT,i,10);
				putdec2(STDOUT,arr[i],10);
				putlin(s3,STDOUT);
				putc(NEWLINE, STDOUT);
				exit(1);
			}
			if(arr[i] >= MPBASE) {
				// negative digit
				putdec2(STDOUT,code,5);
				putdec2(STDOUT,d,10);
				putdec2(STDOUT,i,10);
				putdec2(STDOUT,arr[i],10);
				putlin(s4,STDOUT);
				putc(NEWLINE, STDOUT);
				exit(1);
			}
		}
	}
	if(arr != 0)
		free(arr);
}

/**
 * 	length = MPITOC(mpint,string,size);  convert mp integer  mpint
	to character string  string  of decimal digits of length MPITOC <= size.
 */
int mpitoc(INTEGER n, char** str,int size) {
	char *big 	= " int too big for mpitoc";
	char *s = 0;
	int l = size2(n);
	int rs = 0;
	if(l > MPDIM*10)
		error(big);
	s = mpz_get_str((char *) 0, 10, n);
	rs = strlen(s);
	if(rs > size) {
		s = s + rs - size;
		rs = size;
	}
	*str = s;
	return rs;
}

// convert string of decimal digits at str[i] to mp integer out
// and increment i
void mpctoi(char *str, int i, INTEGER out) {
	// str is a null-terminated string
	while(str[i] == BLANK)
		i++;
	mpz_set_str(out, str+i, 10);
	i+=size2(out);
}

/**
 * MPODIV: Remove largest power (2**L) of 2 from A (so A is both I/O).
 * Function value MPODIV = L.  But MPODIV = -1 if A = 0.
 */
int mpodiv(INTEGER a) {
//	int i,d,j,l,a1,rs;
	int rs;
//	char *e = "error in MPODIV";

	if(compare2(a,0) == 0) {
		return -1;
	}
	// first remove the largest power of MPBASE
//	int *a2 = 0;
//	int s = 0;
//	convert(a,&a2,&s);
//	i = size - 2;
//	a1 = size - 1;
//	if(i>0) {
//		d = a1-i;
//		for(j=2;j<d;j++) {
//			l = j+i;
//			a2[j] = a2[l];
//		}
//		a1 = d;
//		j = 0;
//		d = a2[0];
//		while(d%2 == 0) {
//			j++;
//			d = d/2;
//		}
//		rs = MPBITS*i+j;
//		if(j>0) {
//			if(div1(a,pow(2,j),a) != 0) {
//				mpchk(550,a);
//				error(e);
//			}
//		}
//	}
//	if(arr != 0)
//		free(arr);
	return rs;
}
// C = greatest odd common divisor of A and B
void mpgocd(INTEGER a, INTEGER b, INTEGER c) {
	int i,j;
//	CHAR junkc;
	INTEGER junk1;
	mpcopy(a,c);
	mpcopy(b,junk1);
	i = mpodiv(c);
	if(mpodiv(junk1) < 0) return;
	if(i < 0) {
		mpcopy(junk1,c);
		return;
	}
	i = compare(c,junk1);
	while(i!=0){
		if(i<0) {
			sub(junk1,c,junk1);
			j = mpodiv(junk1);
		} else {
			sub(c,junk1,c);
			j = mpodiv(c);
		}
		i = compare(c,junk1);
	}
}

void mpower(int b, int e, INTEGER m, INTEGER a) {
	INTEGER b2;
	init(b2);
	mpow(b,e,b2);
	mod(a,b2,m);
	clear(b2);
}

int jacobi(int pa, int qa) {
	int ja=1;
	int p,q,r,itmp;

	char *e = "bad argument for JACOBI";

	if(qa <= 0 || pa < 0 || qa % 2 == 0) {
		putdec2(STDOUT, pa,5);
		putdec2(STDOUT, qa,5);
		putc(NEWLINE,STDOUT);
		error(e);
	}

	if (qa == 1 || pa == 1) return ja;

	p = pa;
	q = qa;

	if(q <= p)
		p=p%q;

	if(p == 0) {
		ja=0;
		p=1;
	}

	while(p != 1) {
		itmp = 1 + (p&3);
		switch(itmp) {
		case 1:
			p=p/4;
			break;
		case 2:
			r = q % p;
			if(r == 0) {
				p = 1;
				ja = 0;
			} else {
				q = p;
				p = r;
			}
			break;
		case 3:
			p = p/2;
			r = q%8;
			if(r==3 || r == 5)
				ja = -ja;
			break;
		case 4:
			if(q%4==3)
				ja = -ja;
			r = q%p;
			if(r == 0) {
				p = 1;
				ja = 0;
			} else {
				q = p;
				p = r;
			}
		}
	}

	return ja;
}

// Additional functions
// convert a large integer into a set of integers with number
// of digits <= DIM - 1
void convert(INTEGER n, int** a, int *size) {
	char *str;
	str = mpz_get_str((char *)0, 10, n);
	int l = strlen(str);
	*size = l/DIM + 1;
	*a = (int* )malloc(sizeof(int)*(*size));
	int i=0;
	int start = 0;
	int j=start;
	int value=0;
	for(i=0;i<*size;i++) {
		value = 0;
		for(j=start;j<start+DIM-1 && j < l;j++) {
			value = value*10 + (((int)str[j]) - 48);
		}
		(*a)[i] = value;
	}
}
